1 Preparations

1.1 Load libraries

1.2 Environment

col.os <- "#414a4c"  # Outer Space
col.rp <- "#7851a9"  # Royal Purple
col.cb <- "#b0b7c6"  # Cadet Blue
col.el <- "#ceff1d"  # Electric Lime
col.rm <- "#e3256b"  # Razzmatazz
scale.col.1 <- c()
theme_set(theme_minimal(base_family = "Candara"))
options(dplyr.summarise.inform = TRUE)

2 Load DataSets

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L01-v3_1.html

Land Price Open Dataset

L01_006/_100: Price (JPY/m\(^{2}\))
L01_007: Price Change (%)
L01_022: District Code
L01_023: City
L01_024: Address

admin.sf <- sf::read_sf(
    "input/Administrative_district_Tokyo_2021/N03-21_13_210101.shp"
) %>% rmapshaper::ms_filter_islands(min_area = 1e8)

lp.sf <- sf::read_sf("input/Land-Price_Tokyo_2022/L01-22_13.shp")

# print(st_crs(admin.sf))  # JGD2011
# print(st_crs(lp.sf))  # JGD2000

admin.sf <- admin.sf %>% st_transform("+proj=longlat +ellps=WGS84")
lp.sf <- lp.sf %>% st_transform("+proj=longlat +ellps=WGS84")
admin.code <- readxl::read_excel(
  "input/AdminiBoundary_CD.xlsx", skip = 1
) %>% rename("District" = "行政区域コード")

3 Preprocessing

3.1 Concat Polygons by District

admin.sf <- admin.sf %>% 
  rename(District = N03_007) %>% 
  select(District) %>% 
  group_by(District) %>% 
  summarise(.groups = "drop") %>% 
  left_join(admin.code %>% select(District, City), by = "District")

3.2 Transform and Rename

city.rm <- c(
    "Ooshima", 
    "Niijima", 
    "Kouzushima", 
    "Miyake", 
    "Hachijyou", 
    "Ogasawara"
)

lp.sf <- lp.sf %>% mutate(
    Price = L01_006 / 1e6,
    Price_Change = L01_007, 
) %>% rename(
    District = L01_022,
    # City = L01_023,
    Address = L01_024
) %>% 
    left_join(admin.code %>% select(District, City), by = "District") %>% 
    filter(!City %in% city.rm)

4 Exploratory Data Analysis

4.1 Descriptive Statistics

lp.sf.stat <- lp.sf %>% select(District, City, Price, Price_Change) %>% 
  group_by(District, City) %>% 
  summarise(
    Price.Mean = round(mean(Price, na.rm = TRUE), 3), 
    Price.Median = round(median(Price, na.rm = TRUE), 3), 
    Price.Std = round(sd(Price, na.rm = TRUE), 3), 
    Price_Change.Mean = round(mean(Price_Change, na.rm = TRUE), 3), 
    Price_Change.Median = round(median(Price_Change, na.rm = TRUE), 3), 
    Price_Change.Std = round(sd(Price_Change, na.rm = TRUE), 3), 
    .groups = "drop"
  ) %>% 
  st_set_geometry(NULL) %>% 
  left_join(admin.sf %>% select(District), by = "District") %>% 
  st_set_geometry(value = "geometry")

view.table(lp.sf.stat %>% st_set_geometry(NULL))

4.2 Price and Price_Change Distribution

p1 <- lp.sf %>% select(Price) %>% 
  ggplot() + 
  geom_histogram(mapping = aes(x = Price), binwidth = 0.2, 
                 fill = col.cb, color = "white") + 
  xlim(NA, 20) + 
  labs(x = NULL, y = NULL, title = expression("Price (1e6 JPY)"))

p2 <- lp.sf %>% select(Price_Change) %>% 
  ggplot() + 
  geom_histogram(mapping = aes(x = Price_Change), bins = 100, 
                 fill = col.cb, color = "white") + 
  labs(x = NULL, y = NULL, title = expression("Price Change (%)"))

p1 / p2

# boxplot

4.3 Choropleth Map: Tokyo

p1 <- lp.sf.stat %>% 
    ggplot() + 
    geom_sf(mapping = aes(), data = admin.sf, 
            fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
    geom_sf(mapping = aes(fill = Price.Mean)) + 
    rcartocolor::scale_fill_carto_c(name = "", palette = "Magenta") + 
    geom_sf_text(aes(label = Price.Mean), color = "white", 
                 cex = 2.5, family = "Candara") + 
    labs(x = NULL, y = NULL, title = "Mean Price (1e6 JPY)") + 
    theme(title = element_text(size = 15), 
          legend.position = c(0.05, 0.35), 
          plot.margin = unit(c(0, 0, 20, 0), units = "pt")) + 
    guides(fill = guide_colourbar(
        barwidth = 0.5, barheight = 10, label.position = "right", 
        ticks.linewidth = 3
    ))

p2 <- lp.sf.stat %>% 
    ggplot() + 
    geom_sf(mapping = aes(), data = admin.sf, 
            fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
    geom_sf(mapping = aes(fill = Price.Median)) + 
    rcartocolor::scale_fill_carto_c(name = "", palette = "Magenta") + 
    geom_sf_text(aes(label = Price.Median), color = "white", 
                 cex = 2.5, family = "Candara") + 
    labs(x = NULL, y = NULL, title = "Median Price (1e6 JPY)") + 
    theme(title = element_text(size = 15), 
          legend.position = c(0.05, 0.35), 
          plot.margin = unit(c(0, 0, 0, 0), units = "pt")) + 
    guides(fill = guide_colourbar(
        barwidth = 0.5, barheight = 10, label.position = "right", 
        ticks.linewidth = 3
    ))

p1 / p2

ggsave("fig/Land-price_Tokyo.jpg", dpi = 300, width = 10, height = 10)
p1 <- lp.sf.stat %>% 
    ggplot() + 
    geom_sf(mapping = aes(), data = admin.sf, 
            fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
    geom_sf(mapping = aes(fill = Price_Change.Mean)) + 
    rcartocolor::scale_fill_carto_c(name = "", palette = "Magenta") + 
    geom_sf_text(aes(label = Price_Change.Mean), color = "white", 
                 cex = 2.5, family = "Candara") + 
    labs(x = NULL, y = NULL, title = "Mean Price Change (%)") + 
    theme(title = element_text(size = 15), 
          legend.position = c(0.05, 0.35), 
          plot.margin = unit(c(0, 0, 0, 0), units = "pt")) + 
    guides(fill = guide_colourbar(
        barwidth = 0.5, barheight = 10, label.position = "right", 
        ticks.linewidth = 3
    ))
p1

ggsave("fig/Land-price-change_Tokyo.jpg", dpi = 300, width = 10, height = 5)

4.4 Choropleth Map: Central Tokyo

p1 <- lp.sf.stat %>% filter(District <= 13123) %>% 
    ggplot() + 
  geom_sf(mapping = aes(), data = admin.sf %>% filter(District <= 13123), 
          fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
  geom_sf_text(aes(label = City), color = col.os, 
               cex = 4, family = "Candara") + 
  labs(x = NULL, y = NULL) + 
  theme(plot.margin = unit(c(0, 0, 20, 0), units = "pt"))

p2 <- lp.sf.stat %>% filter(District <= 13123) %>% 
  ggplot() + 
  geom_sf(mapping = aes(), data = admin.sf %>% filter(District <= 13123), 
          fill = NA, color = col.cb, alpha = 0.1, size = 0.1) + 
  geom_sf(mapping = aes(fill = Price.Median), color = "white") + 
  geom_sf_text(aes(label = Price.Median), color = col.os, 
               cex = 4, family = "Candara") + 
  rcartocolor::scale_fill_carto_c(name = "", palette = "Tropic") + 
  labs(x = NULL, y = NULL) + 
  theme(title = element_text(size = 15), 
        legend.position = c(0.05, 0.3), 
        plot.margin = unit(c(0, 0, 0, 0), units = "pt")) + 
  guides(fill = guide_colourbar(
    barwidth = 0.5, barheight = 15, label.position = "right", 
    ticks.linewidth = 3
  ))

p1 / p2 + plot_annotation(
  title = "Median Land Price", 
  subtitle = "(1e6 JPY)", 
  theme = theme(title = element_text(size = 15))
)

ggsave("fig/Land-price_Central-Tokyo.jpg", dpi = 300, width = 10, height = 20)
p1 <- lp.sf.stat %>% filter(District <= 13123) %>% 
  ggplot() + 
  geom_sf(mapping = aes(), data = admin.sf %>% filter(District <= 13123), 
          fill = NA, color = col.cb, alpha = 0.1, size = 0.1) + 
  geom_sf(mapping = aes(fill = Price_Change.Mean), color = "white") + 
  geom_sf_text(aes(label = Price_Change.Mean), color = col.os, cex = 3) + 
  rcartocolor::scale_fill_carto_c(name = "", palette = "Tropic") + 
  labs(x = NULL, y = NULL) + 
  theme(title = element_text(size = 15), 
        legend.position = c(0.05, 0.3), 
        plot.margin = unit(c(0, 0, 0, 0), units = "pt")) + 
  guides(fill = guide_colourbar(
    barwidth = 0.5, barheight = 15, label.position = "right", 
    ticks.linewidth = 3
  ))

p1 + plot_annotation(
  title = "Mean Land Price Change (%)", 
  subtitle = "", 
  theme = theme(title = element_text(size = 15))
)

ggsave("fig/Land-price-change_Central-Tokyo.jpg", 
       dpi = 300, width = 10, height = 10)

5 Spatial Clustering with SKATER

5.1 Standardize Node 4 Features

lp.stat.scale <- lp.sf.stat %>% 
    as.data.frame() %>% 
    select(-District, -City, -Price_Change.Median, -geometry) %>% 
    scale() %>% as.data.frame()

5.2 Compute Neighboorhod list with spdep::poly2nb

Works fine with MULTIPOLYGON, but not with NA.

lp.stat.nb <- spdep::poly2nb(lp.sf.stat %>% select(geometry))
spdep::card(lp.stat.nb)
##  [1] 5 5 6 6 6 5 7 6 5 4 4 7 7 5 6 6 5 5 3 6 4 3 3 6 8 5 5 4 7 4 6 2 7 7 6 4 5 4
## [39] 7 2 4 2 4 4 5 3 4 5 5 4 2

5.3 Compute costs with spdep::nbcosts

lp.stat.costs <- spdep::nbcosts(
    nb = lp.stat.nb, 
    data = lp.stat.scale,
    method = c("euclidean", "maximum", "manhattan", "canberra", 
               "binary", "minkowski", "mahalanobis")
)
# cat(class(lp.stat.costs))

5.4 Making listw with spdep::nb2listw

# unlist(lp.stat.costs)
lp.stat.listw <- spdep::nb2listw(neighbours = lp.stat.nb, 
                                 glist = lp.stat.costs, style = "B")
lp.stat.listw
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 51 
## Number of nonzero links: 250 
## Percentage nonzero weights: 9.611688 
## Average number of links: 4.901961 
## 
## Weights style: B 
## Weights constants summary:
##    n   nn       S0       S1       S2
## B 51 2601 438.3087 2809.367 24129.87

5.5 Find a Minimum Spanning Tree with spdep::mstree and Plot

lp.stat.mst <- spdep::mstree(nbw = lp.stat.listw, ini = 5)
# cat(class(lp.stat.mst))
cat(dim(lp.stat.mst))
## 50 3
lp.stat.mst[1:3, ]
##      [,1] [,2]      [,3]
## [1,]    5   18 0.9166671
## [2,]   18   17 0.7975685
## [3,]   17   19 0.5493890
plot(sf::st_geometry(lp.sf.stat), border = col.os)
# ?plot.mst
plot(lp.stat.mst, coordinates(as(lp.sf.stat, "Spatial")), 
     cex.circles = 0., cex.labels = 0.5, fg = col.os, add = TRUE, 
     col = col.rm)

5.6 Regionalization

5.6.1 ncuts = 2 (ncluster = 3)

lp.stat.reg.2 <- spdep::skater(
    edges = lp.stat.mst[, 1:2], 
    data = lp.stat.scale,
    ncuts = 2
)
summary(lp.stat.reg.2)
##              Length Class  Mode   
## groups       51     -none- numeric
## edges.groups  3     -none- list   
## not.prune     0     -none- NULL   
## candidates    3     -none- numeric
## ssto          1     -none- numeric
## ssw           3     -none- numeric
## crit          2     -none- numeric
## vec.crit     51     -none- numeric
# cat(names(lp.stat.reg.2), sep = "\n")
table(lp.stat.reg.2$groups)
## 
##  1  2  3 
## 29  5 17
p1 <- lp.sf.stat %>% 
    filter(!City %in% c("Oku-Tama", "Hinohara")) %>% 
    ggplot() + 
    geom_sf(mapping = aes(), 
            data = admin.sf %>% filter(!City %in% c("Oku-Tama", "Hinohara")), 
            fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
    geom_sf(mapping = aes(fill = Price.Median)) + 
    rcartocolor::scale_fill_carto_c(name = "", palette = "Magenta") + 
    geom_sf_text(aes(label = Price.Median), color = "white",
                 cex = 3, family = "Candara") +
    labs(x = NULL, y = NULL, title = NULL, 
         caption = "Median Price (1e6 JPY)", size = 1) + 
    theme(title = element_text(size = 15),
          legend.position = c(0.05, 0.35), 
          plot.margin = unit(c(0, 0, 20, 0), units = "pt")) + 
    guides(fill = guide_colourbar(
        barwidth = 0.5, barheight = 15, label.position = "right", 
        ticks.linewidth = 3
    ))

p2 <- lp.sf.stat %>% 
    filter(!City %in% c("Oku-Tama", "Hinohara")) %>% 
    ggplot() + 
    geom_sf(mapping = aes(), 
            data = admin.sf %>% filter(!City %in% c("Oku-Tama", "Hinohara")), 
            fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
    geom_sf(mapping = aes(fill = Price_Change.Mean)) + 
    rcartocolor::scale_fill_carto_c(name = "", palette = "Magenta") + 
    geom_sf_text(aes(label = Price_Change.Mean), color = "white",
                 cex = 3, family = "Candara") +
    labs(x = NULL, y = NULL, title = NULL, 
         caption = "Mean Price Change (%)") + 
    theme(title = element_text(size = 15), 
          legend.position = c(0.05, 0.35), 
          plot.margin = unit(c(0, 0, 20, 0), units = "pt")) + 
    guides(fill = guide_colourbar(
        barwidth = 0.5, barheight = 15, label.position = "right", 
        ticks.linewidth = 3
    ))

p3 <- lp.sf.stat %>% bind_cols(
    tibble(GROUP = factor(lp.stat.reg.2$groups, levels = 1:3))
) %>% ggplot() + 
    geom_sf(mapping = aes(fill = GROUP), 
            size = 0.2, alpha = 0.5) + 
    scale_fill_carto_d(name = "", palette = "Vivid") + 
    geom_sf_text(aes(label = City), color = col.os, 
                 cex = 3, family = "Candara") + 
    labs(x = NULL, y = NULL, title = NULL) + 
    theme(legend.position = "none", 
          plot.margin = unit(c(0, 0, 0, 0), units = "pt"))

p1 / p2 / p3 + 
    plot_annotation(
    title = "Spatial Clustering on Price Stats", 
    subtitle = "", 
    theme = theme(title = element_text(size = 20, family = "Candara"))
    )

ggsave("fig/Spatial-Clustering-Tokyo_n3.jpg", 
       dpi = 300, width = 15, height = 10 * 3)

5.6.2 ncuts = 3 (ncluster = 4)

lp.stat.reg.3 <- spdep::skater(
    edges = lp.stat.mst[, 1:2], 
    data = lp.stat.scale,
    ncuts = 3
)
summary(lp.stat.reg.3)
##              Length Class  Mode   
## groups       51     -none- numeric
## edges.groups  4     -none- list   
## not.prune     0     -none- NULL   
## candidates    4     -none- numeric
## ssto          1     -none- numeric
## ssw           4     -none- numeric
## crit          2     -none- numeric
## vec.crit     51     -none- numeric
table(lp.stat.reg.3$groups)
## 
##  1  2  3  4 
## 29  3 17  2
p1 <- lp.sf.stat %>% 
    filter(!City %in% c("Oku-Tama", "Hinohara")) %>% 
    ggplot() + 
    geom_sf(mapping = aes(), 
            data = admin.sf %>% filter(!City %in% c("Oku-Tama", "Hinohara")), 
            fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
    geom_sf(mapping = aes(fill = Price.Median)) + 
    rcartocolor::scale_fill_carto_c(name = "", palette = "Magenta") + 
    geom_sf_text(aes(label = Price.Median), color = "white",
                 cex = 3, family = "Candara") +
    labs(x = NULL, y = NULL, title = NULL, 
         caption = "Median Price (1e6 JPY)", size = 1) + 
    theme(title = element_text(size = 15),
          legend.position = c(0.05, 0.35), 
          plot.margin = unit(c(0, 0, 20, 0), units = "pt")) + 
    guides(fill = guide_colourbar(
        barwidth = 0.5, barheight = 15, label.position = "right", 
        ticks.linewidth = 3
    ))

p2 <- lp.sf.stat %>% 
    filter(!City %in% c("Oku-Tama", "Hinohara")) %>% 
    ggplot() + 
    geom_sf(mapping = aes(), 
            data = admin.sf %>% filter(!City %in% c("Oku-Tama", "Hinohara")), 
            fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
    geom_sf(mapping = aes(fill = Price_Change.Mean)) + 
    rcartocolor::scale_fill_carto_c(name = "", palette = "Magenta") + 
    geom_sf_text(aes(label = Price_Change.Mean), color = "white",
                 cex = 3, family = "Candara") +
    labs(x = NULL, y = NULL, title = NULL, 
         caption = "Mean Price Change (%)") + 
    theme(title = element_text(size = 15), 
          legend.position = c(0.05, 0.35), 
          plot.margin = unit(c(0, 0, 20, 0), units = "pt")) + 
    guides(fill = guide_colourbar(
        barwidth = 0.5, barheight = 15, label.position = "right", 
        ticks.linewidth = 3
    ))

p3 <- lp.sf.stat %>% bind_cols(
    tibble(GROUP = factor(lp.stat.reg.3$groups, levels = 1:4))
) %>% ggplot() + 
    geom_sf(mapping = aes(fill = GROUP), 
            size = 0.2, alpha = 0.5) + 
    scale_fill_carto_d(name = "", palette = "Vivid") + 
    geom_sf_text(aes(label = City), color = col.os, 
                 cex = 3, family = "Candara") + 
    labs(x = NULL, y = NULL, title = NULL) + 
    theme(legend.position = "none", 
          plot.margin = unit(c(0, 0, 0, 0), units = "pt"))

p1 / p2 / p3 + 
    plot_annotation(
    title = "Spatial Clustering on Price Stats", 
    subtitle = "", 
    theme = theme(title = element_text(size = 20, family = "Candara"))
    )

ggsave("fig/Spatial-Clustering-Tokyo_n4.jpg", 
       dpi = 300, width = 15, height = 10 * 3)

5.6.3 ncuts = 4 (ncluster = 5)

lp.stat.reg.4 <- spdep::skater(
    edges = lp.stat.mst[, 1:2], 
    data = lp.stat.scale,
    ncuts = 4
)
summary(lp.stat.reg.4)
##              Length Class  Mode   
## groups       51     -none- numeric
## edges.groups  5     -none- list   
## not.prune     0     -none- NULL   
## candidates    5     -none- numeric
## ssto          1     -none- numeric
## ssw           5     -none- numeric
## crit          2     -none- numeric
## vec.crit     51     -none- numeric
table(lp.stat.reg.4$groups)
## 
##  1  2  3  4  5 
##  3  3 17  2 26
p1 <- lp.sf.stat %>% 
    filter(!City %in% c("Oku-Tama", "Hinohara")) %>% 
    ggplot() + 
    geom_sf(mapping = aes(), 
            data = admin.sf %>% filter(!City %in% c("Oku-Tama", "Hinohara")), 
            fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
    geom_sf(mapping = aes(fill = Price.Median)) + 
    rcartocolor::scale_fill_carto_c(name = "", palette = "Magenta") + 
    geom_sf_text(aes(label = Price.Median), color = "white",
                 cex = 3, family = "Candara") +
    labs(x = NULL, y = NULL, title = NULL, 
         caption = "Median Price (1e6 JPY)", size = 1) + 
    theme(title = element_text(size = 15),
          legend.position = c(0.05, 0.35), 
          plot.margin = unit(c(0, 0, 20, 0), units = "pt")) + 
    guides(fill = guide_colourbar(
        barwidth = 0.5, barheight = 15, label.position = "right", 
        ticks.linewidth = 3
    ))

p2 <- lp.sf.stat %>% 
    filter(!City %in% c("Oku-Tama", "Hinohara")) %>% 
    ggplot() + 
    geom_sf(mapping = aes(), 
            data = admin.sf %>% filter(!City %in% c("Oku-Tama", "Hinohara")), 
            fill = col.cb, color = col.cb, alpha = 0.1, size = 0.1) + 
    geom_sf(mapping = aes(fill = Price_Change.Mean)) + 
    rcartocolor::scale_fill_carto_c(name = "", palette = "Magenta") + 
    geom_sf_text(aes(label = Price_Change.Mean), color = "white",
                 cex = 3, family = "Candara") +
    labs(x = NULL, y = NULL, title = NULL, 
         caption = "Mean Price Change (%)") + 
    theme(title = element_text(size = 15), 
          legend.position = c(0.05, 0.35), 
          plot.margin = unit(c(0, 0, 20, 0), units = "pt")) + 
    guides(fill = guide_colourbar(
        barwidth = 0.5, barheight = 15, label.position = "right", 
        ticks.linewidth = 3
    ))

p3 <- lp.sf.stat %>% bind_cols(
    tibble(GROUP = factor(lp.stat.reg.4$groups, levels = 1:5))
) %>% ggplot() + 
    geom_sf(mapping = aes(fill = GROUP), 
            size = 0.2, alpha = 0.5) + 
    scale_fill_carto_d(name = "", palette = "Vivid") + 
    geom_sf_text(aes(label = City), color = col.os, 
                 cex = 3, family = "Candara") + 
    labs(x = NULL, y = NULL, title = NULL) + 
    theme(legend.position = "none", 
          plot.margin = unit(c(0, 0, 0, 0), units = "pt"))

p1 / p2 / p3 + 
    plot_annotation(
    title = "Spatial Clustering on Price Stats", 
    subtitle = "", 
    theme = theme(title = element_text(size = 20, family = "Candara"))
    )

ggsave("fig/Spatial-Clustering-Tokyo_n5.jpg", 
       dpi = 300, width = 15, height = 10 * 3)

5.6.4 ncuts = 9 (ncluster = 10)

lp.stat.reg.9 <- spdep::skater(
    edges = lp.stat.mst[, 1:2], 
    data = lp.stat.scale,
    ncuts = 9
)
summary(lp.stat.reg.9)
##              Length Class  Mode   
## groups       51     -none- numeric
## edges.groups 10     -none- list   
## not.prune     0     -none- NULL   
## candidates   10     -none- numeric
## ssto          1     -none- numeric
## ssw          10     -none- numeric
## crit          2     -none- numeric
## vec.crit     51     -none- numeric
table(lp.stat.reg.9$groups)
## 
##  1  2  3  4  5  6  7  8  9 10 
##  3  2 16  1 21  1  4  1  1  1
p1 <- lp.sf.stat %>% bind_cols(
    tibble(GROUP = factor(lp.stat.reg.2$groups, levels = 1:3))
) %>% ggplot() + 
    geom_sf(mapping = aes(fill = GROUP), 
            size = 0.2, alpha = 0.5) + 
    scale_fill_carto_d(name = "", palette = "Vivid") + 
    geom_sf_text(aes(label = City), color = col.os, 
                 cex = 3, family = "Candara") + 
    labs(x = NULL, y = NULL, title = NULL) + 
    theme(legend.position = "none", 
          plot.margin = unit(c(0, 0, 0, 0), units = "pt"))

p2 <- lp.sf.stat %>% bind_cols(
    tibble(GROUP = factor(lp.stat.reg.4$groups, levels = 1:5))
) %>% ggplot() + 
    geom_sf(mapping = aes(fill = GROUP), 
            size = 0.2, alpha = 0.5) + 
    scale_fill_carto_d(name = "", palette = "Vivid") + 
    geom_sf_text(aes(label = City), color = col.os, 
                 cex = 3, family = "Candara") + 
    labs(x = NULL, y = NULL, title = NULL) + 
    theme(legend.position = "none", 
          plot.margin = unit(c(0, 0, 0, 0), units = "pt"))

p3 <- lp.sf.stat %>% bind_cols(
    tibble(GROUP = factor(lp.stat.reg.9$groups, levels = 1:10))
) %>% ggplot() + 
    geom_sf(mapping = aes(fill = GROUP), 
            size = 0.2, alpha = 0.5) + 
    scale_fill_carto_d(name = "", palette = "Vivid") + 
    geom_sf_text(aes(label = City), color = col.os, 
                 cex = 3, family = "Candara") + 
    labs(x = NULL, y = NULL, title = NULL) + 
    theme(legend.position = "none", 
          plot.margin = unit(c(0, 0, 0, 0), units = "pt"))

p1 / p2 / p3 + 
    plot_annotation(
    title = "Spatial Clustering on Price Stats", 
    subtitle = "", 
    theme = theme(title = element_text(size = 20, family = "Candara"))
    )

ggsave("fig/Spatial-Clustering-Tokyo_n2-10.jpg", 
       dpi = 300, width = 15, height = 10 * 3)